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ABSTRACT 

The cell-vertex schemes due to Ni and Jameson, et al., have been 
subjected to a theoretical analysis of their truncation error. The analysis 
confirms the authors^ claims for second-order accuracy on smooth grids, but 
shows that the same accuracy cannot be obtained on arbitrary grids. It is 
shown that the schemes have a unique generalization to axisymraetric flow that 
preserves the second-order accuracy. 
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1, Introduction 


I Most Euler codes written for computing inviscid aerodynamic flows have 

made use of the finite-volume formulation. An early example was described by 

i 

i Rizzi and Inouye [1]. Essentially, they are generalizations of the second- 

j 

order accurate Lax-Wendroff [2] algorithm to nonrectangular mesh geometries. 

I Conservation (which is a necessary ingredient toward ensuring convergence to 

correct weak solutions [3]) comes from regarding the numerical representation 

! 

I of the flow as comprising average states within each computational cell. The 

i 

I rules for updating these values employ interface fluxes which transfer the 

I 

' mass, momentum, and energy between neighboring cells, leaving the total quan- 

1 

I titles present unchanged except for boundary effects. There seems to be 

I little analysis to indicate how accurate these schemes are on nonuniform 

grids. They are clearly second-order accurate on uniform, rectangular grids, 

I 

I and it has been generally assumed that second-order accuracy will still be 

^ obtained on a "sufficiently smooth" irregular grid. 

I Finite-volume schemes based on central-differencing algorithms have 

become very popular following the work of Jameson et al. [4]. Because these 
schemes update the solution via a five-cell rather than nine-cell stencil, 

i 

j their accuracy is easier to analyze and has been treated theoretically by 

i 

I Turkel [5] , whose results have led to modified schemes offering better 

accuracy* Numerical experiments confirming this are reported by Turkel et al. 

[ 6 ]. 

Yet another possibility is the scheme recently introduced by Ni [7]. 
Here the numerical values represent states found at the corners (vertices) of 
the computational cells . Ni originally described his scheme in finite- 
difference terms, with the values approximating point samples of the continuum 
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solution. Subsequently, it has been described in finite-element terras with 
the vertex values implying a bilinear interpolation function within each cell 
(Davis et al. [8], Lohner et al. [9]). This latter view enables a precise 
interpretation of conservation to be given. At a time when various schools of 
thought are learning to use each other's language and ideas, we prefer to 
classify Ni's method under the neutral term of "cell vertex" scheme. One 
attraction of Ni's scheme (which has been very significantly refined by Hall 
[10]) is that it seems to handle arbitrary, even completely unsraooth, meshes 
in a very natural way. However, we will show that claims of uniformly second- 
order accuracy are based on an Incorrect argument. The object of this paper 
is to provide a correct analysis of the accuracy of cell-vertex schemes. 

This is made easier by the fact that these schemes, in a sense, factorize 
the process of updating the solution. Each cell is examined to see if the 
fluxes around it are in balance; if not, some changes are made. The updating 
process at each vertex uses nine vertex states so that, overall, the scheme is 
as complex as Lax-Wendroff . However, the test for balance involves only four 
vertices, making the first stage very amenable to analysis. This first stage 
is the only part to be analyzed here, but this is the only analysis needed to 
determine the accuracy of the steady state which is reached when all cells are 
in equilibrium. That is to say, the integral 

^ (F dy - G dx) (1.1) 

evaluated by the trapezium rule, vanishes around every cell (Figure la). Note 
that there are generally more vertices than cells so that this condition 
leaves the solution undetermined. Boundary conditions need to be added, but 


carefully. If too many are prescribed, then the integral cannot vanish every- 
where. Paisley [11] discusses these isues, including the role of boundary 
conditions in suppressing unwanted "checkerboard” modes in the solution. See 
also Morton and Paisley [12]. 

For problems that are geometrically complicated, or algorithms which 
employ local mesh refinement in difficult areas, it is attractive to use 
meshes composed of triangles, perhaps without any regular structure. Jameson 
et al. [13] have proposed a scheme which also makes use of flow variables 
located at vertices. Again, it is unclear whether to regard the values as 
point samples or as defining a finite-element interpolant (this time 
linear). Clearly, it is not now possible to require that the flux integral 
vanishes around the perimeter of every cell because, in a triangular mesh, 
vertices outnumber cells by roughly 2:1; so the problem would be overcon- 
strained even without boundary conditions. Instead, Jameson^s procedure is 
such that the solution ceases to change when the flux integral vanishes around 
all control volumes such as the one shown in Figure lb. 

The analysis contained in this paper applies equally to Ni*"s and to 
Jaraeson^s schemes because it evaluates the accuracy of the approximate flux 
integral in terms of the geometry of the polygon around which it is 
evaluated. The claim, mentioned earlier, of second-order accuracy arising 
from the trapezium rule is fallacious because, as we refine the grid, the con- 
tour around which we integrate does not remain fixed but shrinks with the grid 
and always contains just a few intervals. In fact, it is easy enough to make 
an analysis based on a local two-dimensional Taylor expansion which demon- 
strates that the integrals are only evaluated to second-order accuracy if the 
polygons meet certain geometric criteria. These criteria are not met by any 
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triangle or by any quadrilateral except parallelograms. They are satisfied by 
all regular n-gons (n > 3) and in a seemingly sporadic manner by various 
irregular n-gons (n > 4). In practice, the errors will be acceptable if the 
polygons have sufficiently small error constants, the formulae for which are 
given. 

There are two offshoots to our analysis. One is that we derive a con- 
sistent finite-element interpretation of the cell-vertex schemes. We show 
that it gives results which differ from those of the finite-difference inter- 
pretation by terms of the same magnitude as the truncation error inherent in 
either. Finally, we consider the application of these schemes to problems 
which are pseudo-two-dimensional, for example, axially symmetric. We show 
that out of a family of plausible extensions, only one member is capable of 
retaining second-order accuracy near the axis. 


2. Description of Cell-Vertex Schemes 

We consider steady inviscid two-dimensional compressible flow governed by 
the Euler equations 


where 


where 



( 2 . 1 ) 


( 2 . 2 ) 


= p[h - i (u^ + v^)] 


P 


Y 


(2.3) 


and h is total specific enthalpy, assumed constant. We intend to satisfy 


these equations in an approximation to their integral form 

f (|1 dy " G dx) = 0 


(2.4) 


where the integral is around any arbitrary closed contour. To this end, the 
plane is divided into cells that are either regularly arranged quadrilaterals 
(Figure la) or unstructured triangles (Figure lb). In either case, the 
vertices are numbered, and a fluid state u^, where ^ = (p, pu, pv), is 
assigned to the ith vertex. Ni's scheme for quadrilaterals approximates (2.4) 
by trapezium rule integration around each cell. This is 


^ell 


I It <1. + I«) - y,) - 7 (a, + G.) ( 


sides 




(2.5) 


where the summation is over all sides of the cells, and all differences are 
taken anticlockwise. Note that D has dimensions of the time derivatives of 
the actual conserved quantities (not their densities), so that £ indicates 
how masses, etc., associated with the points A, B, C, D are to be changed. 
Different time-marching schemes result from specifying different weights with 
which the changes are distributed over A, B, C, D. Ni's original recipe [7] 
reproduced a type of Lax-Wendroff algorithm, and Hall [10] has presented a 
Runge-Kutta version. Here, however, we are interested only in the steady 
state which we assume is characterized by the vanishing of for every 
cell. In practice, this will only happen if the boundary conditions are 
specified carefully. (For example, on a grid of n x n cells there will be 
n^ equations and (n+l)'^ unknowns; hence 2n+l boundary conditions must be 
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prescribed.) Also, the residuals do not in general vanish near a captured 
shock [11]. This implies that our interest here is chiefly in the accuracy 
obtained in smooth parts of the flow. 

In Jameson's triangular grid scheme [13], the average flux across each 
line joining two vertices is evaluated as 

5aB ■ I '"a ^ Ib> ‘i-B - - I <=A * =B> <-B ' ’‘a^’ 

Then masses, etc. , proportional to — added at X and sub- 
tracted at Y. When this has been done for all edges in the mesh, then every 
point has been changed by a net amount proportional to the flux across the 
perimeter of the polygon formed by all triangles meeting at that point. At 
the steady state, we again assume that all these flux integrals vanish - due 
care being taken at the boundary. 

It is the ability to characterize steady state solutions in this rather 
simple way that makes the accuracy of cell-vertex schemes much easier to 
analyze than Lax-Wendroff schemes where the statement of local equilibrium 
typically involves nine fluid states and thirty-two mesh coordinates. The 
statement of equilibrium for a Jameson-type scheme on quadrilaterals involves 
only five fluid states, but analysis of its accuracy needs twenty-four mesh 
coordinates . 


3 . Finite-Difference Analysis 

We wish to know how accurately the formula 


J kA ^ Ib> '»B->’a> 


(3.1) 


approximates the integral 


(3.2) 


Ip = ^ F dy 


taken around the same polygonal contour. At this stage, the values of 
are thought of as point samples, in the spirit of finite difference methods. 
This is most easily done by assuming that ^ is representable by a Taylor 
expansion about some arbitrary nearby origin. We consider each side of the 

polygon in turn. The numerical approximation to the integral is 


I (y, - yA> - ^A> ll * T <»A *B> Ix ^ T <yA * ly 


12 ? 122 
+ y (x. + x„ ) F + 7 - (y . Fn ) P 

4 A B —XX 4 A B — yy 


7 <Va *ByB> Ixyl 


(3.3) 


+ 0 (h") . 


The leading terras in the exact integral are found by introducing a local 
coordinate g varying linearly from 0 to 1 over AB, so that a typical point 
on AB is defined by 


x(5)=CXg + (l-5)x^ 


y (5) = 5 Yg + (1 - 5) y^ 


and 
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F (O = F + X (?) + y (O Fy + y x^(?) 


+ x(?)y(?) F + Y y^O F + O(h^) 


(3.4) 


Inserting these expressions into (3.2) and integrating w.r.t. ? 
yields 


// Fdy - (yj - y^) [F . ‘ ^ (y^ + ^ 


12 2 
+ -T (x. + x.x_ + x_ ) F 
6 A A B B -XX 


* 6 <^Va * ’‘a>’b ’SJ'a * 


xy 


* i <>’a^ * Vb * Jyyl ' 


(3.5) 


Subtracting (3.5) from (3.3) gives the error on AB as 


((Xb - . 2(x,- x^) (yj - y^>4y + (y^ - yf 4^] 


Hence the total error may be arrived at simply by summing the contribu- 
tion from all sides. Clearly we can find the error in the G integral 
similarly. The final result is 


/ F dy = I F Ay - E 2 F^^ 

f G dx ®y GAx— E| G 

) — ^ — 1 — w 


+ 2E, 


F 

-xy 


+ 2E, 


G 

— xy 


E, F 

4 _yy 


E- G 
3 -yy 


(3.6a) 


(3.6b) 
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where 

E, . ^ I (ix)^ 

E 2 = j2 I 

(3.7) 

" l2 ^ Ax (Ay)^ 

= ~ I 

It is pleasing that the error constants take these simple forms and use- 
ful that there are only four of them. Observe that the quantities we wish to 
approximate are O(h^) and that, superficially, the error constants are 
O(h^), making the method first-order accurate. However, some cancellation 
takes place inside the summations in (3,7), which holds out the prospect of 
better accuracy, depending on geometric properties of the grid which we now 
investigate. 


4 , Polygons with Vanishing Error Constants 

It is easy to see that any polygon with central S3nnmetry will have Ej^ 
to identically zero because contributions from opposite sides will 

vanish. Also, given one error-free polygon, we can form others by permuting 
the sides and forming mirror images. It is straightforward to show that an 
error-free polygon remains error-free under rigid rotations. Certain other 


results are also obtainable. 
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Theorem 1 : A quadrilateral Is error-free Iff It is a parallelogram . 

Proof ; For a quadrilateral ABCD, the error constants can be manipulated 
to read 

■ ’^c^ ^^A ■ *B 

Ez = 2S (x^ - Xg + - Xp) + (x^ - x^) (Xg - Xjj) (y^ " Fg ~ 

(4.1) 

“ 2^^^A ■ ^B ^ ^^A ■ ^^A ‘ ’^B ^ ~ 

^4 = ^^A ■ ^^B ■ ^^A ■ ^B ^ 

where 

2S = (x^ - x^) (yg - yp) - (Xg - x.) (y^ - y^,) 

and S Is the area of the quadrilateral. The only terms which may vanish in 

these expressions are (x^-Xg+x^-Xg) , ^y^”yB‘*‘yc~yD) • These quantities, 
divided by two, are the components of the vector which joins the midpoints of 

the two diagonals. Thus to E^ vanish iff the diagonals bisect each other, 

which only happens if the quadrilateral is a parallelogram. 

Corollary : No triangle is error-free. 

Proof: Since triangles are degenerate quadrilaterals, they are excluded 


by Theorem 1 . 
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Theorem 2: Any regular n-gon (n > 3) Is error-free . 

Proof ; For n even the result is trivial by symmetry. A proof, which 
includes n odd, follows by considering the polygon in the complex plane. 
Without loss of generality, one side can be made parallel to the x-axis, and 
each side is then represented by one of the n*-^ roots of unity, say 
(k = 0, ...n-1) where 


2irk 2Trk 

z, = cos + i sin 

k n n 


then 


I 


Y 2nk . 2 2Ttk 2nk . 2irk . 2 2irk . . 3 2trk 

I cos — — + 3i cos — — sin 3 cos sin i sin 

^ n n n n n n 


Y 2iTk ^ 'x ^Tfk 2 2iTkv 

2 cos + 3i cos sin 3 cos (1 - cos ) 

^ n 


n 


n 


. . 27rk .. 2 2iTkv 

- i sin (1 - cos ) 

n n 


= 4(Ej + iE^) 

= (by a similar argument) - 4(E^ + iE^) • 


But 


k=n-l 


I 

k=0 


k=n-l 

I = 

k=0 

zf -1 


3k 


3 
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The numerator vanishes because = !• So the expression equals zero 

unless the denominator also vanishes; that is, unless n = 3 (when the expres- 
sion equals 3). So unless n = 3, E^= E^= 0* 

For n > 4, we can find a variety of asymmetric error-free polygons. The 
general procedure is to prescribe (n-2) vertices. Then requiring that E^ 

to E^ vanish gives four equations for two unknown coordinate pairs. Unfor- 
tunately, these are simultaneous cubic equations for which no systematic 
solution is apparent. To give a flavor of the possibilities, we consider a 
pentagon with prescribed vertices at (0,0), (0,1), and (1,0). We maintain 

symmetry by supposing that an error-free pentagon exists with additional 
vertices at (l+a,b), (b,l+a). Then we require 

^ (Ax)^ = (Ay)^ = 1 - (1+a-b)^ + a^- b^ = 0 

and 

5] (Ax)^Ay = (Ay^)Ax « ab^ - (1+a-b)^ - ba^ = 0. 

Condition (4.2) simplifies to 

3(a-b)(a+l)(b-l) = 0; (4.4) 

if (a-b) = 0, there are no solutions to (4.3). If b = 1, then (4.3) yields 

3^2 ^ 

a + a - a « 0 

a = 0, i (-1 ± /5), 


(4.2) 

(4.3) 
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whereas, if a = -1 then 

3 2 

b - b - b = 0 
b=0, y(l±/5). 

Although there appear to be six solutions, they actually appear in 

pairs. In each pair, the sides are permutations of the same set of vectors. 
The solution to the problem consists of three independent sets of vectors, 
each of which, by permutation, generates 24 different pentagons, if we allow 
the given vectors to permute also. However, out of a given set of vectors, 
only two orderings will produce a convex polygon, that is, those which are 

monotonically ordered by orientation. Those two are mirror images of each 
other. The independent convex pentagons derivable from (4.2) and (4.3) are 

shown in Figure 2. The square is a degenerate case for which one side 

vanishes. 

Error-free polygons can also be constructed by uniting one or more error- 
free polygons that have edges in common; but, though the construction of 
error-free polygons has some mathematical interest, it seems unlikely that any 
practical mesh could be constructed from them, except for completely regular 
meshes. We turn, therefore, to the more practical question of finding 
polygons for which the error constants are acceptably small. 


5. Polygons with Small Error Constants 

The only complete results that have been obtained relate to quadri- 
laterals, for which the error constants are given by (4.1). The local trunca- 
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tion error will be OCh^) provided these expressions are O(h^). We can 
clearly identify factors which are O(h^), and the only possibility for 
second-order accuracy is that 

° (5.1a) 

and 

(7a ■ ~ ° 

Although these conditions are not met by arbitrary quadrilaterals, they 
are met by quadrilaterals which are produced as part of an analytic grid. 
Suppose we have grid generating functions 

X = X (5 , n ) 

y ” y (5, n) 

and the mesh points are the intersections of 5 = mAC , n * nAq , for integer 
values of m, n. Then, if the mesh is systematically refined, A5 , Aq tend 
to zero, and 

2 

^*A " "'C ■ <5. 2a) 

2 

■ >'b >"0 ■ 

provided the cross-derivatives are not large. 

Note that since the error terms Involve only the local intrinsic geometry 
of each cell, neighboring cells may be of markedly different shape, size, or 
orientation. That statement cannot be made for finite-volume schemes 
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[5,11]. Indeed, the local quality of cell-vertex schemes makes them second- 
order accurate on grids generated more straightforwardly than those described 
above. Suppose a coarse grid is determined arbitrarily, and fine grids are 
constructed from it by repeatedly bisecting the cell edges. Within each of 
the original cells, PQRS, there is a mapping defined by 

X = (l-O(i-n) Xp + c(l-n)x^ + Xp + (1-c) nxg (5.3a) 

y = (l-c)d-n) yp + c(l-n)yq + y^ + (l-C) nyg (5.3b) 

which meets the conditions (5.2). 

Related results have been found by Paisley [11], who finds that in the 

2 

cell-vertex method, it is allowable to perturb a rectangular mesh by 0(h ) 

whereas with Jameson's [4] scheme, the allowable perturbation is only 0(h). 


6. A Finite-Element Approach 

As mentioned in the introduction, the schemes under study have been 
described elsewhere as finite element schemes. Whether this is justifiable or 
not, there is one objection that can be removed. This objection is that the 
formula 

f(F dy - G dx) = 1 I (F^ + F^^^) (y^^^ - y^) - (G^ + (x^^^ - x^) 

is obtained by integrating over each triangular (quadrilateral) element 
assuming a linear (bilinear) distribution for both F and G, which is incon- 
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sis tent because there is a nonlinear functional relationship between and 
G. A consistent finite-element model can be derived by assuming that the 
vector ^ varies linearly (bilinearly) over the elements, where 

( 6 . 2 ) 

These variables were introduced by the author to facilitate the con- 
struction of an approximate Riemann solver [14]. They have the property that 
every component of “» ^ bilinear in the components of w^. Thus 



F 


"l”2 


1:1 hw^ 

2y 1 


^ Y+1 2 

+ J IJ 

2y 2 


”2"3 


Y-1 

2y 




(6.3) 


wViere h is a constant for steady flow, equal to the specific total 
enthalpy. In computing a typical term of (6.1) such as 


I 


AB 



there will be subterms of the form 


AB 

ij 


r B 


dy 


The consistent evaluation of this term is 


I 


AB 

ij 


(6.4) 
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(yg - y*) 


2 <^6 ■ >’a> * "iA"jA> 


-6%- >'a> <”iB - "lA> <”JB - "JA>- 


(6.5) 


The first terra contributes the 'naive' approximation (y - y*)(F + F.); 

the second term Is part of a correction that makes the approximation self- 
consistent. The leading contributions to It are due to first derivatives of 
w. The correction term In (6.5) Is, to leading order. 


•r (Yr> ~ Ya)^^- w. + i (y„ - y*)^(x_ - x.) [w. w. + w. w. ] 
6 •'B ■'A ly jy 6 •'B •'A ' B A ^ ly jx lx jy 


* 7 <J'b • >'a> ‘-B - ’'A>Sx“jx' 


When we sura such terms over every side of the polygon, expressions arise 
proportional to the constants Ej , E 2 , E 3 , E^ given by equations (3.7). 
Where these error constants are small enough, the error due to Inconsistency 
Is negligible, and both the consistent and Inconsistent formulae approximate 
the flux Integrals to O(h^). 


7 . An Integral Formulation of Axlsymmetrlc Flow 

Flow which Is Identical In every plane through some axis of symmetry Is 
governed by differential equations that can be written In either of the forms 
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u ^ + G ^ = (v/r) u (7.1) 
or 

(ru)^ + (rF)^ + (r G)^ = p (7.2) 

where u^, F^, ^ are the usual conserved quantities in two-dimensional flow, 
except that r is now a radial coordinate, and v a radial velocity. Here 
p^ is the vector (0, 0, p, O)*^. 

As a basis for numerical work, (7.1) is tempting because it requires only 
slight modification to a two-dimensional code. However, (7.2) is a better 
basis for shock-capturing because it expresses conservation more exactly. To 
see this, imagine that we are performing a truly three-dimensional computa- 
tion, using the same grid in every axial plane, so that each 3-D cell is a 2-D 
cell rotated through some angle 0 with respect to the axis (Figure 3). 
The exact angle does not matter, so long as we avoid the extremes Q and 2 ti. 

If we make a flux balance on such a control volume, the contributions 
from the four faces which do not lie in the radial planes can be evaluated 
analytically, by integrating with respect to the angle. The result is 

“y// / (r£dr-r Gdx) - //jpdA = 0 (7.3) 

where is the interior of the 2-D cell, and its boundary. 

In order to check that a proposed solution was conservative (and hence a 
correct weak solution), we would check that all such integral equations were 
satisfied. A numerical method which imposes this as a constraint will be a 
correct shock-capturing method. In this way, we can give an interpretation to 
conservation of **radial momentum” via the third equation in (7.3). 


I 
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Equatlon (7.3) can be the basis for extending to axisymmetric flow either 
of the cell-vertex methods under discussion. In Ni"s scheme, one integrates 
the quantities r^, r^ around quadrilaterals; in Jameson's scheme, the 
quantity r(^Ay - GAx) is evaluated for each side and used to update the 
adjoining vertices. The problem we now address is how to evaluate these 
things numerically. 


8. Error Analysis for the Axisymmetric Methods 

In either method, the central ingredient is an expression of the form 

I 

I 

/^®r _F dr or G dx. 

I 

i 

I Consider the first of these. Two plausible computational formulae are 


• I * ■'b^' 


( 8 . 1 ) 


^AB 4 ^B^ 


( 8 . 2 ) 


As in the two-dimensional case, we find what the exact value of the 
integral should be in terms of a Taylor expansion about an arbitrary nearby 


point. The result is 


-20- 


^AB " W e^^^A^A ^A^B ■*■ ^B^A ^’^b’^B^-x 


^ T '■'a' ^ ■^a’^b * 'B'>Ir 


■" 24 <^Va •" 4- r^Xg- + TjX^'t 2t3X^X3 f Sr^x^*)!^ 


(8.3) 


■" T2 <^”'a^’‘a ■" 2''a‘'b’‘a ■" "b\ ■" 'b\ * 'a^"b 2'a''b’'b + 3'"b^’‘b> Ixr 


* J <'a^ ■" "a^'b + "a'^B^ + 'b^> irrl * ••• • 


If we make similar expansions for the approximations (8.1), (8.2), and 

subtract from them (8.3), we find error estimates as follows. 


^AB ^AB ■ ^'b'^A* U^'b^'a^^^B *a'-x '*' 6'^'^B " 'a* ^-r 


^(xB-x^)(3r^XA + r^Xg - r^x^ - Sr^x^) 


■ 17<''b'''a> “Va ■ "a’^ + ''b^ ■ ^■'b’'b> ixr 


<8.4) 


<VV'’^B^- Irrl * ••• 


‘aB ■ ^AB ■ ('■b-'a' [- TT<'b"''a><’‘b‘’‘a> Ix ■ T2 <''b‘''a'^ It 


* -n ‘’‘b-*a> <’'a’‘b 


B A —XX 


(8.5) 


■" n <’^b-'a> < Vb 


r„x.) F ] + 
B A — xrJ 
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Corabining terms for a particular polygonal control volume, we obtain the 
following equations relating the exact integral to the approximations (8.1), 
( 8 . 2 ). 

/' r F dr .1/21 4 * 

’ ‘41 (rj^-r/) [F^ + 4 ] - E 34 - E^F^ r ... (8.7) 

wherein we again find the error constants defined in (3.7). At first sight, 

this is a satisfactory result showing that the errors are O(h^) on a smooth 

grid. However, the quantity being approximated is f f rF^ dA, which is 
2 

0(rh ). The problem is that, near the axis, we find cells for which r is 

0(h), and the error terms are 0(h) relative to the true values. On a non- 
smooth mesh, they are even 0(1). 

We observe from (8.6) and (8.7) that these particular errors could be 
eliminated by choosing the linear combination 


AB 


112 
T ^AB ^ T 


AB 


T <'b ■ Ia 'a 4* 'b 4 * ^'b 4'- 


( 8 . 8 ) 


Note that this is the combination obtained if we assume that r, F^, both 
vary linearly along AB, and then integrate exactly. The error analysis for a 
single side now gives 





r^) (r^ + r^) (Xg - x,)‘’F_+ 4^(x„-xJ(r, 


A —XX 


1 

12 ' 


B A' 


^A> 


F 

— xr 


(8.9) 


24 ^’^B 


rX F ] 
A — rr 
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Therefore, the total error around a particular polygon is 


I ^AB " / ^ r(Ax)^Ar + 1 r (Ax) (Ar)^ 


(8.10) 


^Tl^rr J 


r (Ar) + 


The analysis is now identical to the two-dimensional case except for the 
appearance of the mean radii (r) in every term. These factors, however, 
complicate matters considerably. For example, it is no longer true that a 
regular polygon has vanishing error constants because contributions from 
opposite sides no longer cancel in general. To put the error constants into 
forms which make their order of magnitude self-evident needs considerable 
manipulation, which has only been carried out for quadrilaterals. The easiest 
terra to analyze then is 


I - r/ - 2r|r^ . 2r,r/ - r/ 


''c*' ' 'b'' 


* ■ ’^c‘‘ 


+r^-2r^r +2rr^-r^ 

A D AD D 


■ - V - -^B^) - - r^) (r/ - r/) 
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2('a ' '^c)''^D - 'b>[? <'d ■" 'b>^ * 7 ‘■'d - - 7 <’^A + - 7 'V 


3 ^ 
= 2 






’^B ^ 




(’^B 


^’^A '*’ ’^B “ ’^C ” ^’^A 


’^B - ’^C ^ ’^D^' 


( 8 . 11 ) 


The first terra is O(rh^) on grids which are smooth in the sense of 

the previous discussion; it is 0(rh ) on non-smooth grids. The second 

term appears to be O(h^). That it shall be O(rh^) imposes a natural con- 
dition on the grid near the axis. If the cell nearest the axis has one side 
on the axis, then it can be shown that one of the factors ("^A '*’ ^B “ ~ 

rp), (r^ ~ ’^B ~ *^0 '*■ 0(h); the other is 0(rh) — which one is which 

depends on the labelling. It is assumed that this cell is also smooth in the 
same sense as the others. 

We will omit details of the calculations that show the other terms in 
(8.10) to be O(rh^), simply exhibiting the final rearrangements. We find 


I ^(Ax)^Ar = (xg- Xp) (x^ - x^) (r/ - r^^ + r^^ - r^^^) 

+ (x^ - x^) (rg2 - rjj2) (x^ - xg + x^ - x^) 

+ (xg - Xg) (r/ - r^^) (x^ - Xg + x^ -Xg). 


( 8 . 12 ) 
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Here the last two terms are O(rh^) on any smooth grid, and the first 
one is so if the axis is included. Also, with great patience, we can obtain: 


I ^(Ax)(Ar)2 “ Y ’'b ’'C ~ " V 


+ S (r^ + rg + r^ + rj^Kr^ - + r^ - r^) 


^ 2 ^ ^’^A ’^B - - ’^B - 


(8.13) 


^ T ^^^A - ^B ^ ^ ^"b - 


where S is the cell area. 

On smooth grids, the first three lines are O(rh^); the fourth is 

0(h6). 

Now we turn to the terras involving G_ and in equation (8.3). Both 

o 

of these terras are 0(h), but near the axis they alraost cancel so the 
requirement for accuracy is the same as for the other terras. A straightfor- 
ward evaluation of fjp^ dA cannot avoid errors O(h^) so some subtlety is 
needed. 

Integration by parts shows that 

JJp_dA- J r£dx-JJ r • dA , (8.14) 

a Q 

and we can use this to rewrite the terras in G and £ as 

/ r G dx - // £ dA = / r (G - p) dx + // r dA 

8Q a 9a a 


(8.15) 
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where both terms on the RHS are of order (rh^)* Now we can evaluate the 
contour integral by the methods just discussed; the error terms will be 


(G - I r (Ax)^ ^ ■ P^xr^ r(Ax)^ Ar + 

T2 “ P^rr ^ (Ar)^. 


(8.16) 


The summations multiplying 2xr* £rr have already been shown to be 
0(rh^^' first summation can be evaluated as 

I ;(Ax)^ ^ - ^ + + + rjj)(x^ - x^)(Xg - x^Kx^ - x^ + x^ - x^) 

^^A + ^B “ ^’'a - ^B " ''c 

+ ^ S[3(x^ - Xg + x^ - x^)2 + (x^ - Xg - Xg + Xg)2 

+ (x^ + Xg - Xg - Xg)^]. (8.17) 

Every term in this expression can be justified as O(rh^) on a smooth mesh, 
except for the third line. Here S is O(h^) and the first term inside the 
bracket is O(h^), but of the two remaining terms, one will be 0(r^h2) and 
the other O(h^). Thus, the overall order of the expression must be h^. 
However, this does not destroy the second-order accuracy of the scheme near 
the axis, because in that region (£. “ ^)xx 0(r) so that all contribu- 

tions to (8.16) are actually O(rh^). 
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Lastly, we have to evaluate 


// ^ -gY • (8.18) 

a 

A formula for this can be found using finite element techniques, if the 
gradient is rewritten as 


9p _ lr3p 3p 

aTr " J ^ ^ 


where J is the Jacobian of the transformation (C ^ (x, y). Since we 
also have dA = J d? dy, (8.18) can be found straightforwardly. The result 
is 

- [pg(3r^ + 4rg + 3r^ + Zr^^) - p^(3r^ + 2r^ + 3r^ + Ar^^)] 


24 


^Pa ^^B ^ ~ Pc^^B ■" ^ ■ 2^A^^ 


- ^C> 


^ [Pb<^A + " Pd^^A + *C ■ 


(8.19) 


The error term turns out to be 


(r* + - r^)(x* - x„ + 


8 A B 


C - " ^C^^V^D^Pxx 


-27- 


Prr^ 

^ ■" ’^B ■" ’'c ■" ’^D^^^SCx^ - Xg + - Xj^) + 

+ (Xg - Xp) (x^ - x^) (r^ - rg + - r^^)} 

( 8 . 20 ) 

* <'a ■ 'C> <'b ' '^D> '*4 ' *B ■ ’'C ■" *D> ‘*4 * *B ‘ *0 ‘ *D>' '’xr' 

All terras in this expression are O(rh^) provided the grid is smooth and 
contains the axis • The scheme which has emerged, although more cumbersome 

than one would hope for, is the unique cell vertex scheme retaining second- 
order accuracy in axisymmetric flow# It may be noted that it has also one 
property commonly demanded of a numerical code# It is satisfied identically 
by uniform flow parallel to the axis# For such a flow, ^ is a constant, - 

vanishes, and 9p/9r vanishes# Under these circumstances, all contribu- 

tions to the balance equations vanish identically. 

The expressions which have been found as approximations, in effect, to 
3_F/9x and ^G/9r have been derived by other authors from other considera- 

tions. Margolin and Adams [15], in the context of a Lagrangian scheme, sought 
numerical approximations that would make the equation 

^ ~ = div u (8.21) 

an identity at the discrete level# Here V is the volume of a moving parcel 
of fluid. The unique solution to this problem gave the differentiation 

operators derived here. Margolin and Adams report greatly improved accuracy 

in their computations. 
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The same operators were also obtained by Holm et al. [18] as the unique 
operators to preserve in the discrete formulation certain Hamiltonian pro- 
perties of the differential equations. 

If the scheme were to be used in practice, it may be noted that some 
simplifications of the expressions are possible, at least for quadrilateral 
control volumes. The terras involving ^ can be made more compact by 
rearrangement. The terms involving p in (8.19) can be used partly to cancel 
the terms in - ^). What emerges from this is the following, which we 

claim is the simplest version of the unique formula 






Z f2r^ + ’•a£b ‘B ‘ 2^B £b^ 




F <*D " 2c 2 d 2c * 


( 8 . 22 ) 


- T ‘"a • *t>) ^ ^ + 2r^ G^l 




* 51 <Pa ■ £c> ‘V • ('b ■ Pd’ * *C>’ 


^ (Pb - £o> K*a ■ *C> <'b * ''d> - <^A ■ 'c> <*B * “ 
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A striking feature is that the approximations to | r G dx and // £, dA 
are precisely those that would have been obtained by finite-element integra- 
tions assuming bilinear variations for ^ and However, neither of those 

terms individually is second-order accurate. The combination is, given suit- 
able conditions on the grid. 


9. Conclusions 

The cell vertex schemes studied in this paper have been claimed as 
second-order accurate on arbitrary grids. Certainly the results show an 
impressive convergence as the grids are refined [10]. However, we have shown 
that on arbitrary grids, the local truncation error may be first-order. 
Explicit error constants associated with grid geometry have been derived. 
These indicate what sort of cell shapes to avoid and could also be used as 
part of a mesh refinement strategy. 

The extension to axisymmetric flow has been made, and it seems harder to 
achieve accuracy then. The unique formula which retains second-order accuracy 
near the axis has been derived. However, the extension to fully three- 
dimensional schemes is an extremely difficult piece of algebra which has not 
yet been attempted. 
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! APPENDIX 

Some Properties of the Jameson Formula 

[ At convergence, Jameson's scheme for triangles [13] satisfies 

I I ^ ^^i ^i+l^^^i+1 ■ ^i^ " T ^ ^®i ^i+1^ ^*i+l" *i^ “ ° 

I 

around the perimeter of every polygon formed from the set of triangles meeting 
at a common vertex. This has been presented as a set of flux balance 

equations around an overlapping set of control volumes. A straightforward 

I rearrangement of (A.l) is 

I 

■ 7 1 ■ 7 \ ■ Vi> “ “• 

i 

I 

I Now 5] ^y^i+l“ ^i-1^ ” because if the polygon has an odd number of sides, 

j we travel twice round its perimeter to rearrive at our starting vertex, and 

j with even sides there are two closed polygonal paths. Therefore, we can add 

arbitrary constants to the Let us choose ^ which are the 

values at the common vertex (Figure lb). Rewrite (A. 2) as 


I 7 <Io Il> 7 ^^ 1+1 ■ yi-l> - I 7 '5o 7 <’‘ 1+1 - 

The reason for the factors •j is shown in Figure 4, where P, Q are 
the controids of OABS, OBC. Then 


'p " “ T ^^0 ''a ’'b^ ■ T ^^0 "^B “ 7 ^*A ■ ^b’ 


I 
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So (A. 3) can be Interpreted as the flux balance around non-overlapping 
control volumes whose edges join the centroids. This interpretation leads to 
a dual form of the algorithm as follows. 

Scheme II 

For all edges MN in the mesh, compute 

j ^lii In^ ~ T ~ 

where S, T are the other vertices of the triangles adjacent to MN, and use 
this flux to update M, N. 

This has an identical outcome to the standard scheme, which we repeat for 
contrast. 

Scheme I 

For all edges MN in the mesh, compute 

7 - -I <Sm + £n> 

and use this flux to update S, T. 

Indeed, there are two more possibilities, both of which will also yield 
the same results at convergence. 

Scheme III 


For all edges MN in the mesh, compute 
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^ ~ ~ 2^-S 

and use this flux to update M, N. 

Scheme IV 

For all edges MN in the mesh, compute 

j(Fg + F^Xyj^ - y^j) - j(Gg + G^)(x^ - Xjj) 
and use this flux to update S, T* 

Given the interpretation that ^ is piecewise linear within each tri- 
angle, one is tempted to look for control volumes allowing a more accurate 
integration. For example (Figure 5), we may consider the control volume 
formed by joining the centroids of each triangle to the center of each side. 
Along PI the value of F, say, varies linearly from i (F^ + F_) to 

■J ^ 1a ^ with average value ^ ^ ^ view of 

the discussion in Section 6, we neglect the nonlinear variations of ^.) Some 
algebra will show that this does nothing more than reproduce equations (A.l) 
or (A. 3). This result is obvious if we think of integrating over the 

area of the control volumes. Within OAB, is a constant and the area of 

OJPI is one-third that of OAB. 
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